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Abstract 

In  crossing  theory  for  stochastic  processes  the  distribution  of  quantities  such  as  distances 
between  level  crossings,  maximum  height  of  an  excursion  between  level  crossings,  ampli¬ 
tude  and  wavelength,  etc.,  can  only  be  written  in  the  form  of  infinite— dimensional  inte¬ 
grals,  which  are  difficult  to  evaluate  numerically.  A  Slepian  model  is  an  explicit  random 
function  representation  of  the  process  after  a  level  crossing  and  it.  consists  of  one  regression 
term  and  one  residual  process.  The  regression  approximation  of  a  crossing  variable  is 
defined  as  the  corresponding  variable  in  the  regression  term  of  the  Slepian  model,  and  its 
distribution  can  be  evaluated  numerically  as  a  finite— dimensional  integral. 

This  paper  reviews  the  use  and  structure  of  the  Slepian  model  and  the  regression 
method  and  shows  how  they  can  be  used  to  obtain  good  numerical  approximations  to 
various  crossing  variables.  It  gives  a  detailed  account  of  the  regression  method  for  Gaussian 
processes  with  auxilliary  variables  chosen  in  a  recursive  way.  It  also  presents  a  package  of 
computer  programs  for  the  numerical  calculations,  and  gives  numeric.il  examples  on 
excursion  lengths  as  well  as  wavelength  and  amplitude  distributions.  Further  examples  deal 
with  an  engineering  "jump— and— bump"  problem,  and  excursions  for  a  \2— process. 
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Slepian  models  and  regression  approximations 
in  crossing  and  extreme  value  theory 


by 

Georg  Lindgren1  Igor  Rychlik2 

Univ  of  Lund,  Sweden12  and  Colorado  Stale  Univ,  Fort  Collins,  USA2 


1.  INTRODUCTION 

In  the  physical  world,  a  random  function  is  often  described  as  a  sequence  of  local  maxima 
or  minima,  constituting  a  series  of  random  waves.  In  fact,  not  only  the  visual  impression  of 
the  process,  but  also  many  technologically  important  implications  in  such  fields  as  metal 
fatigue  caused  by  random  vibrations,  failure  caused  by  excess  load  on  a  construction,  etc., 
depend  on  the  character  of  the  process  as  a  random  wave  form.  The  basic  objects  in  this 
theory  are  level  crossings  and  local  extremes. 

A  Slepian  model  is  a  random  function  representation  of  the  conditional  behaviour  of 
a  stochastic  process  after  events  defined  by  level  or  curve  crossings.  In  general,  a  Slepian 
model  contains  one  regression  term  with  random  coefficients  which  describe  the  dependence 
on  initial  conditions  such  as  the  slope  at  the  crossing,  the  value  of  the  process  at  a  pre¬ 
determined  point,  etc,  and  one  residual  term,  which  describes  the  deviations  from  the  path 
set  out  by  the  initial  conditions. 

In  crossing  theory  there  are  many  variables  of  great  practical  interest  which  are  very 
difficult  to  analyse,  such  as  the  distance  between  level  crossings,  maximum  height  of  an 
excursion  between  two  level  crossings,  etc.  The  distribution  of  these  quantities  can  usually 
be  written  only  in  the  form  of  infinite— dimensional  integrals  and  there  is  therefore  a  large 
need  for  good  numerical  approximations.  A  regression  approximation  of  a  crossing  variable 
is  defined  as  the  corresponding  variable  in  the  regression  term  of  the  Slepian  model.  Its 
distribution  can  often  be  found  explicitly,  or  expressed  as  a  finite— dimensional  integral, 
which  can  be  evaluated  explicitly. 
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The  model  was  first,  introduced  try  Slepian  (1963)  to  describe  the  behaviour  of  a 
stationary  Gaussian  process  after  a  zero  crossing.  The  term  Slepian  model  was  int  roduced 
by  Lindgren  (1977).  Ditlevsen  (1985a)  gives  a  review  of  different  engineering  applications, 
and  Ditlevsen  (1985b),  (1986),  and  (1988)  contain  further  applications. 

The  regression  approximation  based  on  the  Slepian  model  was  first  used  by  Lind¬ 
gren  k  Rychlik  (1982)  in  connection  with  the  analysis  of  random  waves.  It  has  been  used  to 
find  practically  useful  upper  and  lower  bounds  for  random  wave  characteristics  such  as 
wave  length  and  amplitude  and  related  quantities  often  used  in  reliablility,  ocean  engine¬ 
ering,  structural  and  mechanical  engineering,  and  other  fields  of  technology;  see  e.g. 

Rychlik  (1988,  1989a). 

A  further  important  use  of  the  Slepian  model  is  that  of  asymptotic  expansions.  After 
a  crossing  of  a  very  high  level  many  processes  take  on  an  almost  deterministic  form,  which 
can  be  found  from  the  leading  terms  in  a  Taylor  expansion;  see  Lindgren  (1983,  19S4b, 
1984c),  Ditlevsen  k  Lindgren  (1988). 

Most  of  the  work  on  Slepian  models  and  regression  approximations  in  crossing 
theory  has  been  directed  towards  Gaussian  processes  or  to  functions  of  such  processes. 
However,  the  Slepian  model  is  not  limited  to  Gaussian  processes  even  if  its  structure  may 
not  be  as  explicit  for  non-Gaussian  processes.  The  v^process  is  a  simple  example  of  a 
non— Gaussian  process  for  which  the  Slepian  model  is  analytically  tractable;  see  Aronowich 
k  Adler  (1985),  (1986),  and  Lindgren  (1989). 

In  this  paper  we  shall  give  a  review  of  the  use  and  structure  of  the  Slepian  model  in 
crossing  theory  and  how  it  can  be  used  together  with  the  regression  approximation  to  find 
good  approximations  to  notoriously  difficult  distributions  such  as  wavelength  and  ampli¬ 
tude,  waiting— times  between  crossings  of  constant  and  moving  barriers  and  other  quanti¬ 
ties. 

Section  2  of  this  paper  contains  the  basic  form  of  Slepian  models  after  level  crossings 
and  local  extremes  and  the  simple  regression  approximations  for  the  length  of  an  excursion 
and  for  the  wave  length  and  amplitude  after  a  local  maximum.  In  Section  3  we  give  mathe¬ 
matical  theorems  about  the  interpretation  of  Slepian  models  as  weak  limits  of  conditional 


distributions  given  an  upcrossing  in  the  Kac  L  Slepian  horizontal  window  sense  (see  Kac  & 
Slepian  (1959)),  and  as  empirical  limits  in  t he  observed  series  of  level  crossings. 

Section  4  describes  the  regression  approximation  for  excursion  length  as  well  as  for 
wave  length  and  amplitude,  first  in  a  the  simple  case  with  fixed  regression  variables,  and 
then  in  a  more  complicated  situation  where  the  regression  variables  are  chosen  in  a  recur¬ 
sive  way.  In  Section  5,  finally,  we  apply  the  theory  to  an  engineering  example,  chosen  from 
structural  mechanics,  and  to  the  level  crossing  distances  in  the  \2— process.  Numerical 
examples,  based  on  a  library  of  FORTRAN  subroutines,  are  also  presented  in  this  section. 


2.  INTRODUCTORY  EXAMPLES 


As  first  examples  of  Slepian  models  and  regression  approximations  we  shall  describe  the 
behaviour  of  a  stationary  Gaussian  process  after  uperossings  of  a  fixed  level  and  after  local 
maxima. 

Suppose  £(t)  is  a  stationary  Gaussian  process  with  mean  zero  and  continuously 
differentiable  sample  paths.  When  needed  we  shall  assume  also  the  second  derivative  to 
exist.  Let  the  covariance  function  be  r(t)  =  E(£(s)*£(s+t))  with  spectral  density  R(t), 
i.e. 


and  write 


00 


r(t)  =  f  exp(iurt)R(u;)  du;, 


—oo 


00 


Ao  =  VU(t))  =  r(0)  =  f  RM  du;, 

— 00 


00 


A2  =  V(£'(t))  =  -r"(0)  =  j  uTR(u>)  du;, 

— 00 


Accession  For 

NTIS  GRAM 
DTIC  TAB 

Unannounced 

fin 

w 

□ 

□ 

^  _ 12 

Distribution/ 

_ _J 

Availability 

Codes  | 

Dlst 


iv, 


oo 


A4  =  VU"(t))  =  rlv(0)=  fu>4RMdu>, 


—oo 
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for  the  spectral  moments,  whenever  they  exist. 


Special 


(1) 


For  any  fixed  level  u  we  say  that  a  u-uperossing  occurs  at  time  to  if 
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4(to)  =  u,  4*(t„)  >  n,  (2) 

and  we  write  N{u;  T)  for  the  number  of  u— uperossings  in  the  interval  [0/1  ] •  Then  the 
mean  number  of  uperossings  is  given  by  Rice’s  formula,  (see  Leadbetter  et  al.  1983),  which 
states 

00 

E(N(u;T))  =  T- {  Yf£(0) ^'(0)(u,y)  dy  = 

=  T.fe(0)(u)-E((^(0))+K(0)=u), 

where  z~*”  =  max(0,z).  In  the  normal  case  this  specializes  to 

E(N(u;T))  =  T.(2rr1(A2/A0)1/2exp(-u2/2Ao), 

since  then  £(0)  and  £'(0)  are  independent  and  normal  with  mean  zero  and  variances  Ao 
and  A2,  respectively. 

We  shall  also  need  a  formula  for  the  expected  number  of  marked  uperossings,  i.e. 
crossings  which  satisfy  some  additional  condition,  specified  by  a  set  A.  Let  //( t)  be  a 
random  vector  which  is  jointly  stationary  with  ^(t),  and  write  N(u,A;  T)  for  the  number 
of  u-uperossings  by  £(t),  0<t<T,  which  satisfy  7?(t+-)  e  A.  Then 

E(N(u,A;T))  =  T-/yf.(0)>?,(0)(u,y)-P(^-)6  A  |  ?(0)=u,  i’(0)=y)  dy 


=  T.f?(0)(u)-E(«'(0))+  I{rf(-)eA}  |  <(0)=u);  (3) 

see  Leadbetter  et  al.  (1983). 

We  shall  have  use  for  some  formulas  for  the  conditional  mean  and  covariances  for 
£(t)  given  value  and  derivatives  at  a  specified  point.  Since  Cov(£'(s),£(s+t))  =  -r'(t), 
Cov(£"(s),£(s+t))  =  r"(t),  and  Cov(£"(s),£'(s+t))  =  r'"(t),  one  has 

E(£(s+t)  |  £(s)=u,  <'(s)=y)  =  ur(t)/A0  -  yr'(t)/A2,  (4) 

E(4(s+t)|£(s)=u,  £'(s)=y,  T(s)=z)  = 

=  uA(t)  +  yB(t)  +  zC(t),  (5) 

Cov(£(s+ti),£(s+t2) 1 4(s)=u,  £'(s)=y)  = 

=  r( t2 — tj )  -  r(t|)r(t2)/A0  -  r'(l|)r'(t2)/A2,  ((>) 
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Cov(£(s+ti),£(s+t2)|£(s)=u,  £'(s)=y,  r(s)=z)  = 

=  r( t2 — 1| )  -  r(t|)r(t2)/A0  -  r'(ti)r'(t2)/A2  -  h(ti)h(t->),  (7) 

where 

b(t)  =  Cov(  £(s+t),£"(s)l£(s),£’(s))  =  r"(  t )  -t-(A2/ Ap)r(t) 

“  vu-ouci.eis))1'2  /TT^TT; 


DEFINITION  1:  (Standard  process)  In  the  sequel,  we  shall  illustrate  the  regression 
approximation  by  means  of  a  stationary  Gaussian  low  frequency  white  noise  process  £(t), 
with  covariance  function 

r(t)=^ 

and  spectral  density,  which  is  constant  over  (-/I,  >/3), 

RM  =  1/2V3  for  \u\  <  V3. 

We  shall  call  such  a  process  a  "standard  process".  It  has  V(£(t))  =  V(£'(t))  =  1,  and 

V(Ht))=1.8.  □ 


EXAMPLE  1:  (Slepian  model  after  u-upcrossings.)  Let  the  process  £(t),  t>0,  have  upcros- 
sings  of  the  level  u  at  ti  <  t2  <  ...,  and  consider  the  process  a  fixed  time  t  after  any  one 
of  these  uperossings,  tk,  say.  For  different  tk-points,  £(tk+t)  takes  on  different  values  and 
thus  behaves  like  the  realization  of  a  random  sequence.  The  value  of  £(tk+t)  depends  to 
some  extent  on  £(tk)=u  and  on  (the  unobserved)  £'(tk).  This  dependence  can  be  expressed 
through  the  regression  function,  (cf.  (4)), 

EU(tk+t)  |  £(tk)=u,  £'(tk)=y)  =  ur(t)/A0-yr'(t)/A2, 
while  the  rest  of  the  variation  is  described  by  the  residual  process 


A(t|tk)  =  ((tk+t)  -  E(£(tk+t.)  |  £(tk),  £'(tk)). 
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considering  t  as  a  parameter,  A(t|tk)  becomes  a  realization  of  a  random  process,  which  is 
Gaussian  with  mean  zero  and  a  covariance  function  which,  by  (6),  s  equal  to 

C(t|,  to)  =  r(to-ti)  -  r(t|)r(to)/A0  -  r'(t|)r'(to)/Ao.  (8) 

Thus  we  can  write 


£(tk+t)  —  E(£(tk+t)  |  £(tk),  £'(tk))  +  A(t|tk)  — 

=  4(tk)*r(t)/A0  -  ^(tk)T'(t)/A2  +  A(t|tk).  (9) 


In  formula  (9)  we  could  have  replaced  <f( t k)  by  u,  but  we  have  chosen  to  retain  it 
in  order  to  indicate  the  dependence  on  the  starting  point.  The  derivative  £'(U)  has  still  to 
be  specified,  and  this  is  the  tricky  part  of  the  Slepian  model,  since  £'(ti),  £'(to),  ...  do  not 
form  a  stationary  sequence.  However,  in  the  long  run,  the  empirical  distribution  of  £'(tk), 
taken  over  all  tk,  k=l,2,,..  converges  to  a  Rayleigh  distribution  with  density 


f!"?tk)(y)  =  A^y  exP(-y2/2A2),  y  >  0.  (10) 

Here  ers  stands  for  ergodic,  or  long  run,  distribution,  meaning  that  it  is  the  limit  of  the 
empirical  distribution  of  £'(tk),  i.e.  for  an  ergodic  process, 


#(tk€[0,T|;  C(tl)<Yl 

1  i  rn - 

T-oo  #{tke[0,T]} 


To  describe  the  long  run  properties  of  £(tk+t)  as  tk  runs  through  the  set  of  all 
u— uperossings,  we  therefore  only  have  to  replace  in  (9),  £(tk)  by  u,  and  the  derivative 
£'(tk)  by  a  Rayleigh— distributed  random  variable  ?/  with  density  given  by  (10),  and 
A(t|tk)  by  an  independent  non— stationary  Gaussian  process  A(t)  with  mean  zero  and  co- 
variance  function  C(s,t)  given  by  (8),  thereby  obtaining  a  process 


UO  =  ur(t)/A0  -  /?r'(t)/A2  +  A(t). 


(ID 


This  is  the  Slepian  model  for  £(t)  at  time  t  after  an  uperossing  of  the  level  u. 

One  should  think  of  £„(t)  as  a  stochastic  generator,  which,  by  repeatedly  drawing  reali¬ 
zations  of  t?  and  A(t),  describes  the  long  run  behaviour  of  f(tk+t)  for  k  =  1,  2,  ...  For 
example,  £u'(0)  =  ur'(0)/Ao  -  r)r"(0)/\2  +  A'(0)  =  77  is  distributed  as  £'(tk),  h  =  1,  2,  ... 

□ 


I 
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EX AMPLE  2:  (Slepian  model  after  local  maximum.)  Suppose  £(t),  t>0  lias  twice  contin¬ 
uously  differentiable  sample  paths  and  let  it  have  local  maxima  at  tt  <  t>  <  ...  .  This  is 
similar  to  Exampel  1,  but  now  £'(0  has  its  zero  downcrossings  at  tk,  and  we  want  to 
describe  £(tk+t)  for  different  tk.  One  way  of  doing  this  is  to  let  £(tk+t)  depend  on  the 
height  and  curvature  at  the  maximum  at  tk  through  the  regression 

E(£(tk+t)  |  £(tk)=u,  f(tk)=0,  r(tk)=z)  = 

=  u-A(t)  +  O-B(t)  +  z-C(t). 

Defining  Am(t|tk)  =  £(tk+t)  -  E(£(tk+t) U(tk),£'(tk),£"Uk)),  we  can  write 

£(tk+t)  =  f(tk)-A(t)  +  Htk)-C(t)  +  An(t|tk) 

in  analogy  with  (9).  Both  £(tk)  and  £"(tk)  at  the  maximum  are  random  quantities  with 
a  long  run  distribution  which  turns  out  to  have  the  density  (see  Lindgron,  1970), 

f!(tk),ntk),u,z)  = 

1  9  9 

=  const-(-z)  expf - (Aoz“+2A'>uz+A.4ir)},  z  <  0  .  (12) 

2(AoA4-A3) 

Now,  let  and  (,*m  be  random  variables  with  joint  density  ( 12),  and  let  Am(t) 
be  a  non-stat ionary  Gaussian  process,  independent  of  (£m,  (m),  with  mean  zero  and  with 
covariance  function  Cm(s,t)  given  by  (7).  Thus,  £m,  Cm,  and  Am(  • )  describe  £(tk),  £"(tk), 
and  Am(  •  |tk)  as  tk  runs  through  the  set  of  all  local  maxima. 

The  Slepian  model  for  £(•)  near  local  maxima  is  then  denned  as 

6>ax(t)  =  £mA(t)  +  CmC(t)  +  Am(t).  (13) 

□ 


EXAMPLE  1  (contd.):  (Regression  approximation  for  length  of  an  excursion  above  u.) 

Let  u  >  0  be  a  fixed  level,  and  use  model  (11),  £„(t)  =  ur(t)/A0  —  j] r ' ( t ) / A 2  +  A(t),  to 
describe  the  behaviour  of  £(tk+t)  after  u-uperossings.  Let  T  >  0  be  the  time  of  the  first 
downcrossing  of  the  level  u  for  £u(t),  i.e.  the  length  of  the  excursion  above  u  that 
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started  at  time  0.  It  can  also  be  defined  as  the  smallest  t  >  0  for  which  £u(t)  =  u.  This 
equation  can  be  solved  for  tj,  giving 


V  =  - 


provided  r'(T)  ^  0. 


u(l-r(T)/A0)-A(T) 
r '  (T)/Ao 


However,  this  equation  is  of  no  use  when  we  want  to  find  the  distribution  of  T, 
since  A(t)  is  a  random  process.  Now,  consider  only  the  regression  term  in  £„(t),  which  we 
shall  denote  by  £u(t), 

rf(t)  =  ur(t)/A0  -  t?r'(t)/A2, 

r  »» 
and  define  T  to  be  the  length  of  the  excursion  above  u  by  £„(t),  which  means  that,  i.a., 

£u(Tr)  =  ur(Tr)/A0  -  7?r'(Tr)/ A2  =  u.  (14) 

r 

Again,  solving  for  t),  we  obtain  i]  as  a  function  of  T  only.  Denoting  this  function  by  q, 

n(  l-r(t)  /  A0) 


q(t)  =  - 


r '( t )/ A  2 


we  have 


u(l-r(Tr)/A0l 

1  =  q(Tr)  = - - - . 

r'(T  r )/ A  2 

(Here  we  need  u^O,  since  if  u=0  the  distribution  of  Tr  is  concentrated  at  the  first  zero 
of  r'(t).) 

Now,  there  is  a  simple  relation  between  the  densities  of  ij  and  Tr, 

f  r(t)=f  (q(t))-|J(t)|-I(t),  (1« 

where  I(t)  is  an  indicator  function, 


I(t)  = 


I  if  ur(s)/Ao  -  q(t)r'(s)/A2  >  u  for  all  s,  0<s<t, 
0  otherwise, 


and  J(t)  is  the  Jacobian  of  the  transformation,  i.e. 

dam  (1— r(t)/Ao)r"(t) 

J(t)  —  %  =  u •  { A2/Ao  +  A2- - ^ - }  . 

(M  r'd  r 
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i 


Inserting  the  Rayleigh  density  (10),  we  obtain 

f  r(t)  =  A21(l(t)  exp(-<i(t)2/2A2) •  |J(t)|  -l(t);  (16) 

this  is  called  the  regression  approximation  of  zero  order  for  the  density  of  T. 

To  facilitate  generalizations  we  shall  introduce  the  function 

tf.t(s)  =  |  tf.(t)=u"  =  ur(s)/A0  —  q(t)r'(s)/A2 

obtained  by  replacing  g  by  q(t),  to  make  ^(t)  =  u.  Writing  (14)  as 

G(Tr)/,=  u-ur(Tr)/A0, 

one  can  then  write  the  Jacobian  in  the  general  form 


Figure  1:  Density  Fpr(t)  of  zero  order  regression  approximation  for  length  of  an  excursion 
above  u  for  a  standard  process  £(t);  u  =  0,  1,  2,  3.  i 


FXAM11LF  2  (contd.):  (Regression  approximation  for  wavelength  and  amplitude.)  We 
shall  seek  the  distribution  of  the  wavelength  T  and  amplitude  H  after  a  local  maximum, 
i.e.  the  horizontal  and  vertical  distances  between  a  maximum  and  the  following  minimum. 
Defined  in  terms  of  the  model  function  £max(t)  =  £MA(t)  +  CmC(t)  +  Am(t),  the  wave 
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length  T  >  0  is  the  time  of  the  first  zero  upcrossing  by  ^ax(t),  while  II  =  £max(0)  - 

£max(T). 

The  regression  approximation  of  T,  II  is  defined  as  the  corresponding  quantifies 
Tr,  Hr  in  the  regression  term  £max(t)  =  £mA(t)  +  CmC(t)  of  ^max(t).  Using  A(0)=1, 
C(0)=0,  one  has 

|.A'(Tr)  +  C.O'(Tr)  =  0,  (18) 

Hr  =  U„A(0)  +  („C(0)}  -  (?„A(Tr)  +  <„C(Tr)}  = 

=  f.-(l-A(Tr))-C.-C(Tr),  (19) 

which  can  be  solved  for  To  write  the  solution  in  a  form  that  can  be  generalized  to 

more  complicated  problems,  define 

l-A(t)  C(t)' 

A'(t)  C ' (t) 
and  write  equations  (18)  and  (19)  as 

G(Tr)Um  Cm)T  =  (Hr  0)T 
If  det  G(Tr)  i  0,  this  has  the  general  solution 

(«rn  C.)T  =  G(Tr)-‘  (Hr  0)T  (20) 

Written  explicitly,  this  states  that 
(„  =  Hr-p(Tr)q(Tr), 

C«  =  Hr-q(Tr), 

where 

p(t)  =  -C'(t)/A'(t), 

q(t)  =  -  A'(t)/{( l — A(t))C*(t)  -  A*{ t )C( t ) } . 

r  r 

(The  only  reason  for  writing  the  functions  in  this  form  is  that  the  density  of  T  ,  II  then 
can  be  written  in  a  form  which  is  familiar  in  engineering  literature.) 

The  Jacobian  of  the  transformation  (20)  can  be  expressed  by  means  of  the  function 
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tfh(s)'  =  Mriax(s)'  |  rfax(t)'=0,  ^ax(OKmax(t)=h"  = 

=  hp(t)q(t)A'(s)  +  hq(t)C'(s)  =  h{p(t)q(t)A'(s)  +  q(t)C'(s)} 


p  r 

obtained  by  replacing  £m,  (,'m  by  the  solution  (20),  to  make  T  =t,  H  =h.  Then,  in  analogy 
with  (17), 

^  £th(s)']  t  9 

j(t)  =  -  l-S - 5=i  =  hp’(t)q(t)-. 

(let  CJ(t) 

With  density  of  (m,  given  by  (12),  this  finally  gives  the  following  density  for  Tr,  Hr, 
which  shall  be  called  the  "zero  order  approximation"  of  Lj, 


f  r  r(t,h)  =  f,  -  (hp(t)q(t),  hq(t))- |J(t,h)|  - l(t,h)  = 

'j  1  j|> 

=  const  •I(t,h)*h2|q(t)'V(t)|  • 

2 

•exp  | - ^2  h2q(t)2  (Tm/x)4  {((7r/Tm)2p(t)+l)2  4-  -^jJ,  (21) 

2t 


where  T  =  7r(A2/A4)1/2  is  the  mean  wavelength,  and  e  =  (I-A2/A0A4)1/2  is  a  measure 
of  the  spectral  width  of  the  process.  Further,  I(t,h)  is  the  indicator  function  which  is 
equal  to  one  if  ^tii(s)'  <  U  for  all  s,  0  <  s  <  t,  and  zero  otherwise.  Note  that  if  h  <  0, 
then  £th(s)'  >  0  for  small  s  >  0,  which  implies  that  I(t,h)  =0  for  h  <  0. 


Figure  2:  Level  curves  and  perspectiv  view  of  the  density  f,pr  jjr(t,h)  of  zero  order 
regression  approximation  of  wavelength  and  amplitude  for  standard  process  £(t).  □ 
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The  regression  approximations  described  in  Examples  1  and  2  are  not  very  accurate, 
as  approximations  to  the  true  excursion  and  wave  distributions.  However,  they  can  be  im¬ 
proved  considerably  by  including  some  extra  information  into  the  regression  term,  like  the 
value  of  the  process  or  the  derivative  at  some  predefined  time  after  the  original  crossing. 
This  is  illustrated  in  the  following  example. 


EXAMPLE  3:  (Regression  approximation  of  wavelength  and  amplitude  with  supplemen¬ 
tary  information.)  Consider  the  model  £max(t)  =  £mA(t)  -I-  (mC(t)  +  Am(t)  after  a  local 
maximum,  presented  in  Example  2.  Let  Si  be  a  fixed  timepoint,  write  X]  =  Am  (si),  and 
define  the  new  residual 

A,(t)  =  Am(t)  -  E(Am(t)  I  X,)  =  Am(t)  -  X,b,(t), 

where  bi(t)  is  a  deterministic  function.  Thus, 

£max(t)  =  A(t)  +  CmC(t)  +  Xib|(t)  +  Ai(t),  (22) 

o 

where  Xi  is  normal  with  mean  zero  and  some  variance  o\  =  V(Am'(si)),  and  Ai(t)  is  a 
non-stationary  Gaussian  process,  the  random  variables  (£m,Cm)i  X|,  and  the  process  A((-) 
being  independent. 

We  shall  approximate  the  distribution  of  wavelength  T  and  amplitude  H  by 
simple  functions  of  (£ra,(m)  and  Xj.  To  do  this,  write  H  =  Mq—Mj,  where 

Mq  =  £max(f!)j  M  j  =  ^max(T). 

The  regression  part  of  (22)  is 

£max(t)  =  ^mA(t)  +  (mC(t)  +  Xibi(t). 

r  r 

Define  the  first  order  regression  approximation  Ti  as  the  wavelength  of  £max(t),  i.e.  take 

Ti  to  be  the  first  t  such  that 

fcA'(Tf)  +  C-C'(T[)  +  Xib,'(Ti)  =  0, 

and  define  Mq  =  ^max(O),  Mj  =  £max(Ti),  respectively.  (Note  that  A(0)=1,  C(0)=0, 
b((0)=0,  which  implies  Mq  =  £m.)  Then,  Hr  =  Mq  -  Mj. 
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Thus,  we  have  replaced  the  three  variables  3',  M^,  M(  by  the  approximating  Tr, 
r  r 

Mq,  Mj,  which  are  connected  to  £m,  <;m,  G  through  the  following  system  of  equations, 

Mg  —  i 

m[  =  «.a(t[)  +  <„c(t[)  +  x, inert), 

o  =  ?.A'(T{)  +  CmC'(Ti)  +  X,b,'(Tt). 


Writing 


G(t)  = 


1  0  0 
A(t)  C(t)  b,(t) 
A'(t)  C'(t)  b,*(t) 


M  =  (Mj  Mj  0)T,  X  =  ({„  f.  X,) 1 , 


T 


we  get  the  relations 

M  =  G(Tj)-X,  X  =  G(Ti)-1  -Nl, 


and  finally 


f  r  r  r(mrvnVt)  =  r  v,(u,z,x)- |J(t,mn,m,)| -I(t,h,x),  (23) 

U  1  t)  1 

T  —1  T 

where  (u,z,x)  =  G(t)  •(mQ,m1,0)  are  functions  of  (Gm^nij)  and  hamg-m^ 

In  (23), 

ic  c  X  (u,z,x)  =  ff  r  (u,z)-fY  (x), 

snbSun'M  Snnsm  A1 

where  the  first  density  is  given  by  (12)  and  fv  (x)  is  a  normal  density.  The  Jacobian 

Ai 

J(t,m0,mj)  is  given  by 

^  thx  (  s) '  ]  —  _t 

J(t,mn,m1)  =  -2= - det  G(t)  1  •(uA"(t)+zC"(t)+xb,"(t)), 

u  1  det  G(t) 


where  as  before  4thx(s)'  denotes  the  process  ^max(s)1  with  u,z,x  inserted  to  make  T  =t, 
Mg=mg,  Mj=mj.  The  indicator  function  I(t,h,x)  in  (23)  is  equal  to  one  if  £thx(s)'  <  0 
for  all  s,  0  <  s  <  t,  and  zero  otherwise. 

Finally,  the  marginal  density  of  T  ,H  can  be  obtained  from  (23)  by  numerical 
integration  of  x  over  the  set  { x;  I(t,h,x)  =  1}.  The  accuracy  in  the  approximation 
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depends  on  the  choice  of  Si,  and  there  are  many  possible  suggestions  where  to  put  it;  see 
Section  4  for  further  details,  and  numerical  examples  of  approximations  of  higher  order. 

□ 


3.  SLEPIAN  MODELS  AS  WEAK  LIMITS 

In  Section  2  we  have  presented  two  basic  Slepian  models  for  a  stationary  Gaussian  process 
£(t)  after  level  crossings  and  local  maxima,  respectively, 

Ut)  =  ur(t)/A0  -  J?r'(t)/A2  +  A(t),  (11) 

£max  (t)  —  £mA(t)  +  +  Am(t).  (13) 

We  have  also  given  some  intuitive  interpretations  of  how  they  describe  the  behaviour  of 
4(tk+t)  when  tk  is  an  upcrossing  of  a  fixed  level  u,  or  a  local  maximum,  respectively. 

These  intuitive  interpretations  can  be  made  precise  in  two  different  directions.  One 
is  by  the  horizontal  window  conditioning  principle,  the  other  is  the  ergodic  frequency  intei'- 
pretation.  With  the  horizontal  window  conditioning,  £„(t)  is  defined  as  the  limit  (in  distri¬ 
bution)  as  h  -♦  0  of  £(s+t)  conditioned  on  the  event  that  £(s)  has  a  u-upcrossing  some¬ 
where  in  the  horizontal  window  [tk— h,  tk].  It  was  introduced  by  Kac  &:  Slepian  (1959),  who 
also  showed  the  connection  with  the  frequency  interpretation. 

Both  the  horizontal  window  and  the  ergodic  interpretation  can  be  made  in  the  sense 
of  distributional  convergence  in  a  space  of  sufficiently  smooth  functions.  This  means  that 
also  distributions  of  variables  such  as  excursion  length,  height  of  exceedances,  wave  length 
and  amplitude,  obtained  from  the  Slepian  model,  can  be  interpreted  in  the  horizontal  win¬ 
dow  conditioning  and  in  the  ergodic  frequency  sense,  provided  the  topology  is  chosen  in  an 
appropriate  way.  Such  problems  have  been  dealt  with  by  i.a.  Lamperti  (1965),  Lindgren 
(1977),  Wilson  (1983,  1986,  1988),  Rusakov  &  Seleznjev  (1988),  and  Seleznjev  (19S9a,  b). 
See  also  Leadbetter  et  al.  (1983,  Ch.  10)  for  an  introduction. 

Suppose  £(t)  has  sample  paths  which  are  continuous  functions  (with  probability 
one).  Even  if  not  absolutely  needed  for  most  of  our  results,  we  shall  also  require  the  sample 
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paths  to  be  continuously  differentiable.  A  sufficient  condition  for  this  is  that  the  process  is 
separable  and  the  covariance  function  -r"(t)  of  the  derivative  has  an  expansion 

— r"(t)  =  A2  +  o(| log  |t|  |_0)  as  t  -  0, 


for  some  constant  tv  >  1;  see  Cramer  &  Leadbottcr  (1067).  (In  Remark  3.1  we  comment  on 
what  happens  when  the  derivative  is  assumed  to  exist  only  in  quadratic  mean.) 

Define  (C1,  tf1)  as  the  space  of  continuously  differentiable  functions  on  the  real 
line,  with  metric 


00  _r  S  llP  I  t  I  lX'  (  t)-y'(t)  I 

d(x,  y)  =  |x(0)-y(0)|  +  E  2  r- - lll-r  - , 

r=l  1  +  supj  t|<r|x'(t)-y'(t)| 

>  e-  xn  y  as  n  -♦  oo  means  that  xn(0)  -*  y(0)  and  x^t)  -  y'(t)  uniformly  on  all 
bounded  intervals.  Further,  let  (C»  tft‘)  denote  the  subspace  of  (C»,  tf1),  consisting  of 
functions  such  that  x(0)=u,  and  write,  for  r>0,  Br  =  [0, r] ,  and  B  =  [-r,0].  We 
consider  £(t)  as  a  process  defined  on  the  space  (C«,  tf1),  with  £(t;  x)  =  x(t);  e.g.  £(tk+-) 


is,  with  probability  one,  and  element  of  C C  We  also  use  the  same  notation  for  functionals 
defined  on  (C1,  &1)  and  the  corresponding  random  variables. 

Let,  for  any  interval  B, 


N(u;  B)  =  #{tk  €  B}, 

be  the  number  of  u-uperossings,  that  fall  in  B.  Similarly,  for  any  set  A  6  let 
N(u,A;  B)  =  #{tk  6  B;  ^(tk+-)  6  A}, 

be  the  number  of  u-upcrossings  in  B  for  which  the  sample  function  starting  at  tk 
satisfies  the  condition  given  by  the  set  A.  These  definitions  generalize  the  definitions  of 
Nfu;  T)  and  N(u,A;  T)  in  Section  2. 

We  shall  give  some  examples  of  relevant  sets  A.  The  first  example  descriltes  the 
length  of  an  excursion.  Let  rk  be  the  length  of  the  excursion  above  u  starting  at  tk,  i.e. 

rk  =  inf { t  >  0;  £(tk+t)  =  u}, 

and,  for  x  €  C^,  let 
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t  =  inf { t  >  0;  x(t)  =  u}. 

Then,  for  =  (x  6  CM;  r  <  T}, 

N(u,At;  B)  =  #{tk  e  B;  7-k  <  T} 

is  the  number  of  excursions,  starting  in  B,  with  length  less  than  T. 

As  a  second  example  we  consider  the  height  of  an  excursion  above  u.  With  M  = 

{x  6  CM;  sup  x(t)  <  m}, 

0<t<r 

N(u,M  ;  B)  =  #{tk  €  B;  sup  £(tk+t)  <  m} 

0<t<rk 

is  the  number  of  excursions  in  B  which  do  not  exceed  m. 

We  can  now  define  two  different  distributions  on  the  space  (CM,  tf1).  First,  let  P 

1  v  u’  ir  u 

be  the  Palm  distribution  for  £(  • )  near  u— uperossings,  defined  by 
Pu(A)  =  E(N(u,A;  B,))/E(N(u;  B,)). 

(Note  that  here,  by  stationarity,  Bi  can  be  replaced  by  any  interval  of  length  one.)  Second, 
define  P  T  to  be  the  empirical  distribution  of  £(•)  near  u-uperossings,  i.e.  (for  T  large 
enough  to  make  N(u;  B^)  >  0), 

Pu,T(A)  =  N(u’A;  BT^N(u;  BT}’ 

i.e.  Pu  j  gives  equal  weight  1/n  to  each  of  the  n  =  N(u;  B.p)  functions  £(tk+0 
starting  at  ti,  ..,  tn. 

THEOREM  3.1:  The  Slepian  model  £„(t)  =  ur(t)/Ao  —  /^r ' ( t ) / A o  +  A(t)  has  distribution 
given  by  Pu  considered  as  an  element  in  (C1,  tfl). 

Proof:  The  finite— dimensional  distributions  of  £u(t)  are  given  by  those  of  P„.  This 
follows  from  the  general  formula  (3)  for  the  expected  number  of  marked  crossings.  One  has 
only  to  take  Tj(t)  =  (£(t+St),...,  £(t+sn)).  Since  the  finite-dimensional  distributions  deter¬ 
mine  Pu,  and  £u  has  continuously  differentiable  sample  paths,  the  t  heorem  follows.  0 
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THEOREM  3.2:  If  £  is  an  ergot! ic  process,  then,  with  probability  one,  P  T  converges 
weakly  to  Pu  in  (C^,  tS’p  as  T  -*  oo. 

Proof:  The  theorem  states  that,  with  probability  one,  the  empirical  probability  measure 
over  the  realizations  of  £(tk+t),  k=l,2,...,  in  (C^,  tf^)  converges  weakly  towards  Py  as 
T  -» oo.  It  follows  from  the  the  ergodic  theorem,  (see  Cramer  &  Leadbetter,  1967),  that,  for 
each  A€tfJ,  with  probability  one, 

T— 1  •  N(u,A;  nT)  -*  E(N(u,A;  B|)), 

T—1  •  N(u;  Bt)-E(N(u;  B,)), 

which  shows  that  Pu  j(A)  -»  PU(A),  with  probability  one.  Since  ( C ^ ,  tf>)  is  separable 
with  the  metric  d,  it  follows  from  Billingsley  (1968),  that,  with  probability  one,  Py  r^(A)  -> 
Pu(A)  for  every  A  €  tfC  □ 


THEOREM  3.3:  The  conditional  distribution  of  {£(•)  |  N(u;  B _ j1)>0}  converges  weakly 

to  Pu  in  (Ci,  tfj)  as  h  [  0. 


Proof:  The  finite  dimensional  conditional  distributions  of  {^(si),...,^(sm)  |  N(u;B_^)  >  0} 
converge  to  Pu  as  h  |  0,  i.e. 

Ptt(Sj)  <  yj,  j=l,...,m  |  N(u;  B_h)  >  0}  = 

P(N(u;B_^)>0  and  ^(sj )  <yj,  j=l,...,m} 

—  - - - - -» 

P(N(u;B_h)>0) 


E(N(u,A  ;  B0) 

- - 1 - =  P(A  ), 

E(N(  u  ;  B  i ))  y 

where 


(24) 


Ay  =  {*  €  C(i;  x(sj)  <  yj,  j=l,...,m}. 

To  prove  (24),  let  to  be  the  last  u-uperossing  of  £(t)  before  0,  so  that  to  >  -h  if  and 
only  if  N(u;B_ji)  >  0.  Then  replace  the  event 
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{ N( u; S _ h)  >  °  and  £(sj)  <  yj,  j=l,...,m} 

by 

{t0  > -h  and  £(t0+Sj)  <  yj, 

The  error  will  be  of  small  order.  Then  (24)  is  a  consequence  of  the  general  result,  for 
stationary  point  processes  v ,  that  h-1P(/>(Bh)  >  0)  -*  E(i/(Bt))  as  h  -  0,  applied  to  the 
point  processes  of  u-upcrossings  {tk}  and  to  the  thinned  point  processes  of  u-upcrossings, 
which  satisfy  the  extra  condition  £(tk+Sj)  <  yj,  j=l,...,m. 

To  show  the  full  weak  convergence  on  (0^,  ^)  we  need  the  tightness  of  the 
conditional  distributions  of  £(0)  and  of  the  continuity  modulus  of  £', 

snp  U'(»W(t)| 

| s-t  |  <  6 
s,teB 

given  that  N(u;B_h)  >  0.  In  fact,  P(|£(0)|  >  K  |  N(u;B_jj)  >  0)  -*  0  as  h  -»  0  for  all 
K  >  u,  so  that  £(0)  is  tight. 

Further  we  have  to  show  the  conditional  tightness  of  (£',£),  i.e.  that,  for  each 
e  >  0,  r  =  1,  2, 

lim  1  im  P{cjr  (£',6)  >  e  |  N(u;B  l)  >  0}  =  0. 

<^0  h-*0  _n 

One  can  use  the  same  technique  as  above,  considering  the  point  process  N  ^  of  u- 
upcrossings  tk  such  that 

"Br+1(U)<?''<)  -  e' 

where  Br+i(tk)  is  the  sphere  of  radius  r+1  and  center  tk-  Then 

P{^B  U’,5)  >  e  |  N(u;B_h)  >  0}  <  P{N^(B_h)  >  0  |  N(u;B_h)  >  0} 

which  tends  to  E(N^(B1))/E(N(u;B1))  as  h  -*  0.  Here,  N^(Bj)  <  N(u;Bj)  and 
N^(Bj)  -*  0,  with  probability  one,  as  6  -*  0,  which  proves  the  tightness.  □ 
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Remark  3.1 :  If  £  has  continuous  sample  paths  with  mean  square  derivatives,  one  can 
define  the  model  £u  as  an  element  in  the  space  (0,,,#,,)  of  continuous  functions  (with 
x(0)  =  u),  with  metric  for  uniform  convergence  on  hounded  intervals  and  the  corresponding 
statements  to  Theorems  3.1  and  3.2  hold.  However,  variables  like  excursion  length  and 
height  after  a  u— upcrossing  are  examples  of  variables  which  are  continuous  on  the  space 
(Ci,  tfji)  but  not  on  the  simpler  space  (Cu,  tfu)-  □ 


4.  REGRESSION  APPROXIMATIONS 


4.1  Regression  decomposition 

Assume  £(  • ),  T)(  • )  is  a  pair  of  smooth  processes  defined  on  the  same  probability  space, 
and  let  T  be  the  first  zero  crossing  of  Suppose  that  we  are  interested  in  the  joint 
density  of  the  crossing  time  T  for  £  and  the  value  of  rj  at  T,  i.e. 


(T,H)  =  (T({),  ,(T)), 


(25) 


which  are  sometimes  called  marked  exit  time  and  mark,  respectively.  In  this  section,  we 
shall  illustrate  the  general  regression  approximation  method,  by  constructing  a  sequence  of 
approximations  to  the  joint  density  of  T  and  II. 

The  regression  method  is  based  on  a  decomposition  of  the  processes  £(t)  ,77ft)  into 
linear  regression  terms  and  independent  residual  processes.  More  precisely,  let  (  = 

( Cf  n>l’  be  a  vector  of  regressors  from  the  linear  space  spanned  by  £(  • ),  r?( • )  and 

let  the  joint  density  of  (p...,  Cn  be  a  smooth  function.  Suppose  that  the  processes  £  and 
t)  can  be  decomposed  as  follows, 

n 


£(s)  =  bQ(s)  +#E  Cjbjfs)  +  Aj(s), 


(26) 


n 


r?(s)  =  cQ(s)  + 


E 

i=l 


Cjcj(s)  +  A2(s), 


where  the  functions  bj  and  Cj  are  defined  explicitly,  and  Aj,  A.>  are  zero  mean 
Gaussian  residual  processes,  independent  of  (Cp...,(n)-  Furthermore,  we  shall  assume  that 
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Var((j)  =  1,  for  all  i,  which  only  a  matter  of  scaling.  The  conditional  expectations  of  the 
( £(s ) , tj( s ) )-processes ,  given  C  =  (Cj ,-,Cn)»  are  §iven  by 

E(e(s)|C)  =  bQ(s)  +  £  (^(s)  =  £f(s), 

1-1  (27) 

E(?7(s)  1 0  =  Cn(s)  +  E  C;C:(s)  =  7?r(s). 
u  i  =  l 

In  Example  1  in  Section  2  we  had  only  one  process, 

£„(t)  -  u  =  ur(t)/A0  -  u  -  t?r'(t)/A2  +  A(t), 

1 1° 

which  is  of  the  form  (26)  with  n  =  1,  Ci  =  ??/ V ( 77)  1  “  and 
h0(t)  =  ur(t)/A0  -  n, 
b1(t)=-r'(t)/A2-V(,)1/2, 

A, ft)  =  Aft). 

Example  2  involves  two  processes, 

^u'(t)  =  *.A'(t)  +  CmC'(t)  +  A'(t), 

£u(0  =  CmA(t)  +  CmC(t)  +  A(t), 

and  is  of  the  form  (26)  with  n  =  2,  (j  =  £m/V(£m)1//2,  (2  =  Cm/V(Cm)^“,  and 
b0(t)  =  cQ(t)  =  0, 

h,(t)  =  A'(t).V(U1/2,  l),(D  -  (’'(t.)-V(CJ1/2, 

Cj(t)  =  Ad).V(^)1/2,  c2(t)  =  C(t)-V(Cm)1/2, 

Aj(t)  =  A'(t),  A2(t)  =  A(t). 

We  shall  now  study  the  general  model  (26)  in  more  detail  in  order  to  approximate 

the  distribution  of  wave-length  and  amplitude,  i.e.  T,  H.  The  regression  approximation  of 

r  r 

(T,H)  is  defined  as  the  corresponding  variables  in  the  regression  functions  (£  ,7/  ),  i.e. 

they  are  obtained  by  replacing  in  (26)  the  processes  (£,77)  by  the  regression  curves.  We 

r  r 

shall  denote  the  regression  approximation  of  (T,H)  by  (T  ,11  ),  i.e. 
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(Tr,Hr)  =  (TUr),  r;r(Tr)).  (28) 

As  has  hern  shown  in  the  introductory  examples,  the  general  structure  of  Slepian 
models  for  the  Gaussian  processes  is  the  decomposition  (26).  Although  most  of  our  appli¬ 
cations  of  the  regression  method  is  in  connection  with  Slepian  models,  it  is  not  limited  to 
the  Slepian  models  and  Gaussian  residuals.  However,  since  we  require  that  all  finite-di¬ 
mensional  distributions  of  the  residual  processes  Ap  A9  are  given  in  explicit  form,  the 
residual  processes  are  usually  asumed  to  be  Gaussian  or  functions  of  Gaussian  vector 

processes;  see  Rychlik  (1987d,  1989b).  In  addition,  in  Section  5.2,  we  shall  present  Slepian 
o 

models  for  the  x~~Process-  that  case,  Slepian  models  will  be  defined  by  more  compli¬ 
cated  decompositions  (functions  of  the  linear  decomposition  (26)).  The  regression  method 
can  still  be  used  for  approximation  of  the  (T,H)-distributions  and  it  is  defined  as  the 
corresponding  variables  in  the  nonlinear  regressions  obtained  by  replacing  the  residual 
processes  by  zero. 

Finally,  the  results  presented  in  this  section  can  be  easily  extended  to  more  compli¬ 
cated  marked  crossing  problems,  when  £  and  77  in  (26)  are  vector  valued  processes,  i.e. 

T  =  (Tlv..,Tk)  is  a  vector  of  the  first  k  zero  crossings  of  £;  see  Rychlik  (1987b,  1988). 


4.2  Simple  regression  approximation 

We  now  turn  to  the  evaluation  of  the  density  of  (T1 ,  llr).  Obviously,  T1  and  ll1  are 
functions  of  (Cp---i(n)  alone.  However,  since  we  are  interested  in  the  density  of  (Tr,  Hr), 
we  need  an  inverse  mapping  that  expresses  ((. ,  (.>)  as  a  function  of  (Tr,  Hr,  ), 

leaving  ((j,...,(n)  unaffected.  This  can  be  done,  using  the  implicit  oefinition  of  Tr,  i.e. 

b0(Tr)  +  E  C.bj(Tr)  =  0.  (29) 

u  i=l 

The  required  variable  transformation  can  be  written  as  follows, 


_  00  _ 


C,b,(Tr)  +  C,MTr)  =  -bn(Tr)  -  E  Cil-i(Tr), 

C,c,(Tr)  +  (,c2('l'r)  =  Ilr  -c0(Tr)  -  £  Cici(Tr). 

i — 3 

r 

The  linear  equation  system  (30)  has  a  unique  solution  ((p  (0)  for  T  =t  if  the 
determinant  of  the  following  matrix  is  not  equal  to  zero, 


G(t)  = 


VO  b2(t)' 


[Cj  (t)  c2(t)J 

As  in  Examples  1  and  2,  we  can  then  solve  (30)  and  obtain  (t^,  £.,)  as  a  function  of 

r  r 

(T  ,H  ,  (3,,..,^).  Denoting  this  function  by  (qp  q0),  i.e. 


q,(t,h,x3,...,xn)l  ,  f  -b0( t )-x3b3( t . .-xnbn( t ) ‘ 

=  tl(t) 


Qo(  . .  ,x^) 


■h-c0(t)-x3c3(t)-...-xncn(t)j 


r  I* 

the  conditional  density  of  (T  ,11  |  C  •••,(„)  can  be  obtained  by  a  simple  variable 
transformation,  given  in  Lemma  4.1  below. 

Before  we  state  the  lemma,  we  shall  introduce  two  important  functions,  £ti,x(s)  and 
*?thx(s),  defined  by  "£r(s)  |  £r(t)=0,  r?r(t)=h,  C3,...,<n=x"  and  "r/r( s)|£r(t)=0,  //r(t)=h, 
C3v,Cn=x",  respectively,  i.e.  for  a  fixed  (t,h,x)  =  (t,h,X3,...,xn), 

r  11 

^thx(s)  =  bQ(s)  +  q1(t,h,x)b1(s)+q2(t,h,x)b2(s)  +  E  XjbjCs) 

(33) 

r  n 

i?thx(s)  =  cQ(s)  +  q1(t,h,x)c1(s)+q2(t,h,x)c2(s)  +S  XjC^s). 

For  any  continuously  differentiable  function  f,  let  I(t;f)  be  an  indicator  function, 
defined  equal  to  one  if  f(s)  has  no  zero  crossings  in  the  interval  (0,t)  and  zero  otherwise, 


I(t;f)  = 


1  if  f  ( s)  >  0  or  f(s)  <  0  for  all  s,  0<s<t, 
0  otherwise . 
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LEMMA  4.1:  Let  (Tr,  Hr)  be  the  regression  approximation  of  ('I',  11),  defined  by  (28).  If 

r  p 

the  determinant  det  G(t)  ^  0,  the  conditional  density  of  T  .II  |  0ti---.Cn  (t,h|x)  = 

(t,h|x3,...,x„)  is  given  by 


f(t,h  |  x)  =  I(t;^thx)  *  |  J(t,h,x)  |  f(,^2|(3...(n((l1(t>h,x),  q2(t,h,x)|x), 


(35) 


where  ^  is  the  conditional  density  of  (t,£>  |  (3,...,(n  and  the  functions  qi,  qo 

are  defined  by  (32).  The  Jacobian  is  given  by 


J(t,h,x)  =  - 


-j^£thx  (  S) 


s-t 


det  G(t) 

p  p 

where  £thx(s)  is  defined  by  (33)  and  the  indicator  function  I(t;£thx)  is  given  by  (34). 


(36) 


□ 


p  p 

The  marginal  density  of  (T  ,  H  )  can  be  obtained  by  multiplying  the  conditional 
density  (35)  by  the  density  of  C3,...,Cn  and  integrating  the  x-variables,  i.e.  by  computing 
E(f(t,h|C3 . Cn)). 

The  Jacobian  J(t,h,x)  is  written  in  the  form  (36)  to  indicate  the  explicit 

p  p  p 

dependence  of  the  T  ,11  -density  on  the  derivative  £ti,x(t)\  i.e  on  the  derivative  of  the 

p  p 

regression  £  at  the  zero  crossing  at  time  t,  when  the  mark  7/  takes  the  value  h. 

r  r 

The  density  of  T  ,11  is  related  to  the  marked  zero  crossing  intensity  through  the 
following  alternative  formula: 


fTrHr(t,h)  =  E(f(t,hj<3.-Cn))  = 

=  E(I(ti5r)Ur(t)'|  |  Cr(t)=0,  I|r(t)=h)fcr(l)i,jr(t)(0,l,).  (37) 

J» 

where  I(t;£  )  is  defined  by  (34).  The  proof  of  formula  (37)  is  based  on  the  observation 
that  the  conditional  density  (£r(t),r?r(t)|  C3---Cn)  is  given  by 


f?'(t).1r(t)l<3.-<n(0'h|><)=  |(letG(l)l  lfC.,C2l(3-Cn(ql(t’h'X)’q2(t’h’X)|x)- 
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4.3  An  explicit  formula  for  clensiiv  of  (TJjJ 


It  is  instructive  to  write  the  true  density  of  T,II  in  a  form  similar  to  formula  (35),  based 
on  the  regression  approximation.  To  obtain  this  form,  just  introduce  the  residual  processes 
A1  and  A0  into  the  functions  q(,  q9  and  the  Jacobian  J,  (36),  defined  in  Lemma  4.1. 
More  precisely,  define  new  q— functions  q,  and  a  Jacobian  T  as  follows, 


q^t.hfXg,. 


■xn’A> 


“Vt^”x:Wl^””“xii,>n(t^l(t)' 

h-c0(t)-x3c3(t)-...-xncn(t)-A2(t), 


(38) 


J(t,h,x,Aj)  =  - 


det  G(t) 

where  A(s,t)  =  (A^s),  A9(t))  and  the  function  £u»x(s)  is  defined  by 


(39) 


£thx(s)  =  l>0(s)  +  q1(t,h,x,A)b1(s)+q2(t,h,x,A)b.2(s) 


+  S  x.b.(s)  +  A,(s). 
i=3  1 


(-10) 


Then  the  true  conditional  density  of  (T,H  |  )  can  be  written  as  an  expec- 

tation  of  the  density  (35)  over  all  sample  paths  of  A  =(Aj,A9),  with  the  modified  q—  and 
T— functions,  i.e. 


fT,H|  C3,...,Cr/t,h,X^  ~  EA1,A.,pt;^lhx)‘  l^(l>,1»x^i)l  • 

* f  C»  C2 1 C3  ...Cn^i^t,h,x,A)’  n2M»x»A)lx) 


(41) 


where  I(t;  * )  is  defined  by  (34). 

The  density  of  (T,H)  can  be  written  in  the  same  form  as  (37)  of  (Tr,Hr),  by 
replacing  the  regressions  £r,  if  by  the  original  processes  £,  //,  i.e. 

fT|H(t,h)  =  E(I(t;0 1  ?(t)  |  |  *(t)=0,  i7(t)=h)f^(t)  J?(t)(0,h).  (42) 

Observe  that  formula  (42)  holds  for  a  broader  class  of  processes,  which  contains  more  gene¬ 
ral  nonlinear  decompositions;  see  Durbin  (1987),  Rychlik  (1987c).  However,  for  the  purpose 


-25- 


of  approximation  and  numerical  calculation,  the  explicit  formula  (11),  based  on  the 
decomposition  (26),  is  more  useful  than  formula  (42). 

The  expectation  (41)  is  difficult  to  evaluate  exactly,  since  the  indicator  I(t;£u,x) 
depends  on  the  whole  realization  of  <fthx,  i.e.  on  Aj(-)-  Nevertheless,  the  formula  is  the 
basis  for  improvement  of  the  accuracy  of  the  regression  approximations  of  the  (T,II)- 
density,  by  approximation  of  the  residual  processes  Ap  A.-,,  by  simpler  processes  Ap  A0, 
which  are  functions  of  a  finite  number  of  random  variables  Y  =  (Yp...,Y.  ),  say,  defined 
on  the  same  probability  space  as  Ap  A0.  The  ( ApA9)-processes  can  be  chosen  in  many 
different  ways,  e.g.  as  cosine  polynomials,  splines  etc.  In  the  next  subsection  we  shall  use 
regression  curves. 


4.4  More  comnlica ted  approximations 

Let  Y  =  (Yp...,Y^)  be  a  random  vector  defined  on  the  residual  processes  Ap  A>?.  The 
general  formula  for  the  regression  approximation  (t£>h£)  of  the  (T,  H)-density  is 
obtained  by  replacing,  in  (41),  the  Aj  -  and  A9-residual  processes  by  their  conditional 
expectations,  given  Y,  i.e.  by  A  =  (ApA,2), 

A1(Y)(s)  =  E(A1(s)|Y),  A2(Y)(s)  =  E(A2(s)|Y).  (43) 

More  precisely,  the  regressions  ApA9  in  (43)  define  a  new  decomposition  of  the  (£,77)- 
processes  by  splitting  the  residuals  ApA9,  i.e. 

«s)  =  ?r(s)  +  A1(Y)(s)  +  A1(Y)(s), 

tfs)  =  7,r(s)  +  A2(Y)(s)  +  A,(Y)(s), 

where,  for  i=l,2, 

Ai(Yr..Yk)(s)  =  A.(s)  -  E(A.(s)|Yp..Yk). 

r  r 

The  regression  approximations  T^,  are  defined  as 

(Tk'  Hk'  =  (T(?r(''  +  4,(Y)(-)),  7,r(Tp  +  A,(Y)(Tp). 


(«) 
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The  conditional  density  flv(t,li|x,y)  of  ('1^,11  |  (J,V)  can  he  obtained  in  a  similar 
way  as  in  Lemma  4.1.  The  explicit  formula  can  he  written  as  follows.  With 


qj(t,h,x,y)l  _1  f  -b0(t)-x;3h3(t)-. 

=  G(t) 

q2(t,h,x,y)J  lh— c0(t)— x3c3(t)— . 


•_xnbn(t)_Ai(y)(t) 

•~xncn(t)~A2(y)(t)- 


J(t,h,x,y)  =  - 


^?hxy(s)3s=t 
(let  d(t) 


and  ^thxy(s)  and  7?Ihxy(s)  defined  by 


^thxy(s)  =  b0(s)  +  q1(t,h,x,y)b1(s)+q2(t,h,x,y)b2(s)  (47a) 

+  ^  x.b.(s)  +  A.(y)(s), 
i=3  1  1  1 

^thxy(s)  =  cQ(s)  +  q1(t,h,x,y)c1(s)+q2(t,h,x,y)c.2(s)  (47b) 

+  S  x.c.(s)  +  A.,(y)(s), 
i=3 

the  conditional  density  fk(t,h|x,y)  of  (T|^,H^|(,Y)  is  equal  to 
fk(t,h|x,y)  =  I(t;£thxy) *  |J(t,h,x,y)| 


fClC2  |  C3 ...Cn^qi^t,h,X,y^  Wt,h’X,y^X^  (4S) 

(Observe  that  the  Y— variables  are  independent  of  £.) 

The  density  of  the  simplest  regression  approximation,  formula  (35),  corresponds  to 
the  choice  k=0,  i.e.  there  are  no  extra  regressors  at  all  besides  Thus  the  approximating 
variables  will  be  denoted  as  Tq,Hq  and  the  conditional  density  f(t,h|x),  formula  (35), 
will  be  denoted  by  f^(t,h  |  x). 

The  regressors  Yj,...,Y^  can  be  choosen  with  great  freedom  to  contain  simple  but 
efficient  information  about  (T,H).  One  of  the  main  features  of  the  structure  of  approxi¬ 
mation  (48)  is  that  it  can  easily  be  improved  by  just  adding  more  regressors  in  a  recursive 
way.  In  the  following  subsection  we  present  a  recursive  procedure  of  choosing  the  regressors 
Yi,...,Yk,  suitable  for  inplementation  in  a  recursive  programming  language. 
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4.5  A  recursive  procedure  for  choosing  the  rocressors  Y, . V, 


i  rv 

The  regressors  Yp...,Yk  should  he  chosen  in  such  a  way  that,  the  conditional  density  of 
(Tk,Hk  |0  is  a  good  approximation  to  the  exact  formula  for  the  density  (41)  of  (T,  H), 
with  as  few  regressors  as  possible.  Consequently,  strategies  of  choosing  the  Yj-variables  are 
usually  quite  complicated  and  can  depend  both  on  the  properties  of  the  Ohx~,  and  7/thx- 
functions  and  on  the  residual  processes  A  =  (ApA0).  We  shall  now  present  a  class  of  re- 
gression  approximations  (Tj.,  Hk)  for  the  (T,II)-distribution,  based  on  recursive  selection 
of  regressors  Yj. 

Assume  that  we  have  designed  a  strategy  to  select  a  unique  variable  Y^,  defined  as 
a  function  of  the  residuals,  A(s,t)  =  (Aj(s),A0(t)).  More  precisely,  assume  that  we  have 
been  given  a  zero-mean  Gaussian  process  A,  with  a  covariance  r^  =  (rA,,rA0,rA  A^’ 
and  deterministic  functions  f  =  (fj,f9),  corresponding  to  the  functions  £ti,x,  T) thx,  and 
assume  further,  that  we  have  defined  a  functional,  Proc,  on  the  sample  paths  of  A  =  u;, 
with  f  as  the  parameter,  such  that  the  regressor  Yj  is  defined  by  Yj(u;)  =  Proc(u> ;  f). 

The  same  procedure,  Proc,  can  then  be  used  to  define  further  regressors  Y0, 
provided  it  has  the  following  structure: 

(a)  Y^ProcfAJ)  is  a  Gaussian  random  variable, 

(b)  the  regression  curves  A(Yj)(s,t)  in  (43)  can  be  evaluated  explicitly, 

(c)  for  any  real  y ^ ,  the  residual  processes  A1(y])(s),  A.,(yj)(s)  in  (44)  are 

zero-mean  jointly  Gaussian  processes,  with  explicitly  defined  covariance 

function  rAiV 

A  |  Y,=y( 

Then,  the  regressors  Yj,  ...,Yk  can  be  defined  recursively: 

0)  =  Proc(A;  £thx>  ^thx)> 

(ii)  for  any  value  y  =  (yp.-My^j)  of  the  regressors  Yj,...,Yk_1, 

=  Pr0c( A(y);  £thxy?  f?lhxy)> 

where  A(y)(s)  =  ( A,(s)-E( A, (s)  |  Y, ...Yk_l ),  A.,(s)-E( A2(s) |  Y, ... Yk_, ))  are  the 
residual  processes  (44)  and  £u,xy,  fjihxy  are  defined  by  (47). 
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We  turn  now  to  the  problem  of  choosing  the  procedure  Proc(  •  )■  Many  different  al¬ 
gorithms  can  be  proposed  and  we  shall  classify  them  into  two  main  categories.  The  first 
one,  which  will  be  called  deterministic  procedures ,  selects  Yj  from  the  linear  space  span¬ 
ned  by  A,  independently  of  the  properties  of  £thx,  rfihx-  Consequently,  the  conditional  co- 
variances  yk_y  do  not  depend  on  y.  As  examples,  we  give  the  simple  procedure 

ProCj(A;  £ihxT/thx): 

Choose  Yj  =  Aj(Sj ),  where  Sp  0<s^<t,  is  a  fixed  point 
such  that  Var(Aj(Sj))=  supn<s<t  Var(A](s)), 

and  the  more  complicated  algorithm 

Proc0(A;  £thxi?/thx)-' 

Choose  Yj=Aj(Sj),  where  Sj,  0<Sj<t,  is  a  fixed  point 
such  that  s»Po<s<t  Var(Aj(s)|  AjfSj))  is  minimized. 

The  second  class  of  procedures  for  choosing  the  first  regressor  Yp  will  be  called 

stochastic  procedures ,  and  here  r^  |  Yk=y  ^ePen(^s  explicitly  on  y.  This  happens 

r  r 

when  one  allows  Yj  to  depend  on  £thxi  VtUx-  An  example  of  a  stochastic  procedure,  used 
in  Rychlik  (1987d,  1989b),  is  the  following  algorithm: 

ProCj(A;  ^thx^thx)’ 

Choose  Y^  =  A^(s^),  where  Sp  0<Sj<t,  is  a  fixed  point 
such  that,  if  £thx(0)>0,  P(£thx(Sj)  >  0)  is  minimized, 
otherwise,  i.e.  if  £u,x(0)<0,  P(^i,x(Sj)  <  0)  is 
minimized. 

The  difference  between  a  stochastic  and  a  determinstic  procedure  is  that  in  the 
stochastic  procedure  the  regressors  Yp...,Yjp  depend  on  the  sample  path  of  the  residual 
processes  Ap  A9,  leading  to  more  accurate  approximations. 

r  r 

Finally,  we  give  a  recursive  formula  for  the  conditional  density  (T^H^  |  ) 

defined  by  recursive  selection  of  the  regressors  Y  =  (Yp...,Y^).  Obviously,  using  the 
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k  r  r 

conditional  density  f  (t,h  |  x,y),  given  by  (48),  of  (Tj.,IIj.  |  (,Y),  we  obtain  the 
(T£,h£|  ^-density  in  the  form  of  a  k— dimensional  integral, 

fk(t,h|x)  =  Ev(fk(t,h|x,Y)),  (49) 

which  has  to  he  evaluated  numerically.  Since  the  regressors  Y  are  choosen  in  a  recursive 
way,  the  value  of  the  conditional  density  f^(t,h|x)  depends  only  on  the  functional  Proc, 

and  the  starting  values  A,  Ohx,  ^thx-  We  shall  express  this  dependence  explicitly  by  the 

k  k 

notation  fpr0C(A-£[h  ^(t,h|x).  The  recursive  formula  for  f  (t,h|x)  is  then  obtained 

by  writing  the  integration  in  (49)  in  iterated  form,  i.e. 


=  (  ey2 ...Yt«V.''l*.y,Y2...Yk))  fY,(y)  dy 


=  I  fProc(  A I  Y,=yi{Ehxy,>)tl,«y)(t’h  1  x)  fY,(y)  dy' 

The  main  advantage  of  the  recursive  representation  of  the  conditional  density  of 

j»  j* 

(Tjt,H^|  0  is  that  it  can  be  easily  implemented  in  a  computer  program,  in  particular  when 
the  programming  language  allows  recursive  functions. 


4.6  Rounds  for  the  T.H— density 

In  some  applications,  one  needs  exact  bounds  of  the  approximations  error.  In  the  following 
we  shall  briefly  present  upper  and  lower  bounds  for  the  (T,H)— density,  based  on  regression 
approximations. 

Let  (t,h,x)  be  a  fixed  vector  for  which  bounds  for  the  conditional  density  of 
T,H|  (g,...,(n  are  required.  Suppose  that  Aj(t)  =  A.-^t)  =  Aj(t)  =  0,  a.s.  This  seems  at 
first  to  be  a  strong  assumption,  but  can  easily  be  attained  by  inclusion  of  A^t),  A0(t), 
A|(t)  in  the  vector  ((p...,Cn)-  (Observe  that  this  was  not  required  when  the  regression 
approximations  of  T,II  were  defined.)  Consequently,  formula  (11)  for  the  density  of 
T,H|  (g,...,(n  simplifies  to 
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l(t;£ttix) 


I J(t,h,x) 


fT,H|C3-Cn(t’h|X)  “  EAj 

•^,0!C3...Cn((M(l'h’x)’(12(t’h-x)lx)- 


(50) 


Assume  that  from  knowledge  of  the  values  (t,h,x),  we  can  determine  whether  the 
first  crossing  of  £thx,  defined  by  (40),  is  an  up-  or  a  down-crossing,  i.e  we  assume  that 
(flOU'(O))  is  included  in  (^,...,(n.  Suppose  that  the  first  zerocrossing  of  £thx  is  an 
upcrossing,  and  hence  the  expectation  in  (50)  is  the  following  probability, 

EA|{I(t;&h*)}  =  pUthx(s)  <  0  for  all  s,  0<s<t).  (51) 

Observe  that  £thx(s)  is  a  continuously  differentiable  Gaussian  process.  In  Rychlik  (1987d, 
1989b)  we  have  proposed  a  recursive  procedure,  similar  to  the  regression  approximation,  to 
bound  probabilities  of  type  (51). 


4.6  Numerical  examples 

In  this  section  we  shall  illustrate  higher  order  approximations  to  the  excursion  length 
distribution  in  Example  1  and  to  the  wavelength  and  amplitude  density  of  Examples  2 
and  3.  Computer  programs  for  the  algorithms  will  be  presented  in  the  next  section,  where 
also  details  about  the  approximations  can  be  found. 

First  consider  the  excursion  length  distribution  for  a  standard  process.  Figure  3,a-d, 
shows  excursion  length  distributions  after  upcrossings  of  levels  u  =  0,  1,  2,  3,  calculated 
with  k  =  1,  2,  3  extra  regressors.  For  comparison,  the  figures  also  show  the  simple 
approximations  with  k  =  0.  As  seen,  it  is  necessary  to  use  extra  information  in  order  to 
catch  the  bimodal  distribution  after  a  zero  upcrossing.  For  higher  levels  even  the  simple 
approximation  gives  reasonable  results,  and  in  fact  for  u  =  2,  3  the  curves  for  k  =  1,  2,  3 
coincide. 
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Figure  3:  Regression  approximations  according  to  formula  (10)  to  excursion  length  density 
after  upcrossing  of  a  level  u  for  a  standard  process;  (a)  u=0,  (b)  u=l,  (c)  u=2,  (d)  u=3. 
The  order  of  approximation  is  k  =  0,  I,  2,  3. 


As  a  second  example  we  show  the  higher  order  approximations  to  the  wavelength 
and  amplitude  distributions  considered  in  Examples  2  and  3.  As  can  be  seen  from  Figure  4, 
(a)-{c)  higher  order  approximations  are  smoother,  in  particular  at  the  left,  smaller  end  of 
the  distribution.  The  zero  order  approximation,  shown  in  Figure  2,  has  a  cut  off  point  at 
t  «  4,  to  the  left  of  which  it  is  zero. 
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4.7  Computer  programs 

The  numerical  examples  presented  in  this  paper  have  been  evaluated  by  means  of  a  library 
of  computer  programs,  CROSSREG,  which  contains  the  regression  approximations  for  the 
following  densities: 

EXCREG:  excursion  time  above  the  level  u  for  the  Slepian  model  ^  (t);  Example  1; 

WAMPREG:  wavelength  and  amplitude  after  local  maximum;  Example  2; 

JUMPREG:  length  of  jump  for  a  car  travelling  on  a  random  road;  Example  5.1; 

o  o 

CHIREG:  excursion  time  above  the  level  u“  for  the  Slepian  model  for  the  \  “-process; 

Section  5.2,  formula  (62),  n  =  2. 

The  programs  are  implemented  in  FORTRAN  77  and  have  been  run  on  a  PC/AT  and  also 
on  a  MicroVAXII. 

Using  the  library  one  can  evaluate  both  the  simple  regression  approximations, 
without  any  additional  regressors,  Section  4.2,  and  the  more  complicated  approximations 
discussed  in  Section  4.4.  The  additional  regressors  Y  =  (Yp...,Yj.)  are  selected  recur¬ 
sively,  using  a  sligthly  modified  stochastic  procedure  of  the  type  Proc^;  see  Section  4.5. 

The  main  modification  lies  in  the  introduction  of  a  stopping  criterion,  which  allows  us  to 
use  different  numbers  of  regressors  for  different  sample  paths  of  the  residual  process.  The 
stopping  criterion  tests  whether,  for  given  values  of  t  and  y  —  (y |  the  conditional 

density  of  the  regression  approximation,  frpr  j  y(t|y)  is  sufficiently  close  to  the  conditional 
density  fip|y(t|y)  of  the  true  T.  This  stopping  rule  reduces  the  necessary  amount  of 
numerical  computations  drastically. 

The  programs  can  be  used  on  a  PC/AT,  with  a  reasonable  computation  time,  for  as 
many  as  four  additional  regressors.  As  can  be  seen  in  the  examples  in  this  paper,  this  is 
enough  to  compute  very  accurate  approximations  for  a  large  class  of  processes  of  practical 
interest.  In  addition,  in  order  to  speed  up  computations,  we  have  designed  a  numerical 
integration  algorithm,  specially  well  suited  for  the  regression  approximation.  A  more 
complete  description  of  the  programs  will  be  given  elsewhere. 
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5.  APPLICATIONS 

5.1  Random  vibration  of  mechanical  structures 

Stochastic  processes  and  linear  filtering  theory  is  often  used  to  describe  the  behaviour  of 
mechanical  structures,  such  as  bridges,  high  towers,  aeroplane  wings,  car  suspensions,  etc, 
subjected  to  external  random  forces.  Loads  created  by  the  environment,  such  as  wind  and 
wave  forces,  are  described  as  stochastic  processes,  while  the  structure  itself  is  described  in 
terms  of  linear  or  non-linear  transfer  functions.  Processes  of  interest  are  forces  and  dis¬ 
placements  within  the  structure,  and  the  crossing  and  extremal  properties  of  these  are  of 
great  importance  for  the  reliability  of  the  structure.  Examples  of  important  quantities  are 
high  peak  values,  length  of  excursions  above  high  levels,  and  period  and  wave  analysis. 

Most  models  used  to  describe  random  vibration  are  high  dimensional,  and  the  exter¬ 
nal  forces  even  form  a  continuous  random  field,  as  for  example  the  wind  over  the  whole 
extension  of  a  long  bridge.  We  shall  here  present  a  simplified  example  with  only  one  forcing 
process  and  a  transfer  function  with  a  simple  non-linearity,  which  contains  many  of  the 
characteristic  features  of  a  random  mechanical  system. 

Example  5.1  ("Jumps  and  bumps";  Lindgren  (1981))  Consider  the  movements  of  a  car 
travelling  with  constant  speed  on  a  rough  road,  considered  as  a  locally  stationary  stochastic 
process.  The  suspension  system  of  the  car  should  guarantee  safe  roadholding  and  high  ride 
comfort.  The  roadholding  depends  on  the  forces  between  wheel  and  road  and  the  ride 
comfort  mainly  on  the  acceleration  forces  on  the  driver  and  passengers.  An  undesirable 
event  occurs  when  the  force  between  a  wheel  and  ground  is  less  than  a  required  level.  In  an 
extreme  case,  a  wheel  can  even  leave  the  ground  for  a  short  while  -  it  jumps  and  bumps. 
We  shall  describe  the  jump-and-bump  situation  by  means  of  a  Slepian  model  and  use  the 
regression  technique  to  derive  the  distribution  of  the  length  of  the  flight  starting  at  the 
jump. 

To  simplify  the  analysis  we  shall  consider  a  one— wheeled  car  with  mass  m  which 
drives  with  constant  speed  v  on  a  randomly  profiled  road.  The  suspension  system  consists 


of  a  linear  spring  with  length  P  ,  stiffness  k,  and  damping  coefficient  c.  (A  four-wheeled 
can  be  treated  similarly  as  a  multivariate  system;  see  Jogreus  ( 1 9S6) . )  Let  Y(t)  denote 
the  height  of  the  road  above  a  zero  reference  level  at  time  t,  and  let  X ( t )  be  the  extension 
of  the  spring  from  its  unloaded  position.  (Then  Y(t)  =  Y0(vt)  where  Y0(s)  denotes  the 
road  elevation  at  distance  s  from  the  starting  point.) 

If  the  wheel  has  no  mass,  the  following  differential  equation  governs  the  movements 
of  the  car  as  long  as  no  jumps  occur, 

m(X"(t)  +  Y"(t))  +  cX'(t)  +  kX(t)  =  -gm  (52) 

where  g  =  9.81  m/s2  is  the  acceleration  of  gravity.  This  implies  that  E(X(t))  =  -gm/k. 
As  long  as  no  jumps  occur,  the  vertical  acceleration  of  the  car  is 

X"(t)  +  Y"(t)  =  -  g  -  m_1(cX'(t)  +  kX(t)), 

but  if  the  wheel  is  not  permanently  attached  to  the  ground,  this  acceleration  can  never  be 
less  than  -g.  Therefore  a  jump  occurs  every  time  the  normal  force 

N(t)  =  cX'(t)  +  kX(t)  (53) 

(which  has  mean  -gm)  has  an  upcrossing  of  the  zero  level.  Then  an  excursion  starts, 
during  which  the  movements  of  the  car  and  spring  are  purely  deterministic. 

If  tk  denotes  the  time  of  a  jump,  the  spring  expands  exponentially  during  the  jump 
according  to  the  equation 

X(tk+t)  =  X(tk)  exp(-kt/c),  (51) 

while  the  car  follows  the  parabolic  curve 

PQ  +  X(tk)  +  Y(tk)  +  t(X'(tk)+Y'(tk))  -  gt“/2.  (55) 

The  road  level  Y(tk+l)  is  of  course  not  affected  by  the  jump,  but  the  fact  t  hat  a 
jump  did  occur  at  tk  affects  its  distribution,  so  we  need  the  conditional  distribution  of 
Y(tk+t)  given  that  a  jump  has  occurred. 

We  now  assume,  as  is  often  done  in  technical  literature,  that  the  road  profile  Y(t) 
is  a  Gaussian  process,  which  is  stationary  at  least  locally.  This  means  that,  if  no  jumps 
were  allowed  to  occur  and  the  wheel  was  in  permanent  contact  with  the  ground,  also  the 


-  36  - 


normal  force  N(t)  and  the  spring  extension  X(t)  would  he  Gaussian  processes.  One  could 
then  use  the  theory  in  previous  chapters  to  describe  the  crossings  and  extremal  properties 
of  N(t)  and  X(t).  However,  if  jumps  may  occur,  the  Gaussian  character  of  the  process  is 
destroyed,  but  if  jumps  are  rare,  one  could  think  of  the  process  as  being  renormalized  after 
each  jump,  and  use  the  Slepian  model  based  on  normal  theory  to  describe  the  local  be¬ 
haviour  near  jumps.  The  assumptions  that  jumps  occur  unfrequently  is  a  realistic  assump¬ 
tion  for  real  cars. 

In  this  example,  we  shall  assume  that  jumps  do  not  occur  at  all  and  that,  X ( t )  is 
always  defined  by  equation  (52),  and  we  nail  consider  equations  (54)  and  (55)  as  fictitious 
excursions,  to  be  described  by  Slepian  models. 

We  assume  that  Y(t)  =  Yu(vt)  is  a  stationary  Gaussian  process  with  mean  zero 
and  (speed  dependent)  covariance  function  ry  possessing  spectral  density 

RY(u;)  =  v_1Ryju;/v), 

where  Rv  is  the  spectral  density  for  the  road.  We  let  X(t)  be  the  stationary  (Gaussian) 
*0 

solution  of  equation  (52),  rewritten  as 

mX"(t)  +  cX'(t)  +  k(X(t)+gm/k)  =  -  mY"(t),  (52') 

and  write  N(t)  =  cX'(t)+kX(t).  All  covariance  properties  of  the  involved  processes  can 
then  be  derived  from  Ry(o>)  and  the  transfer  function  11(a))  of  the  system  (52'), 

HM  = - (56) 

k  +  i  c  u>  —  moT 

We  need  the  following  variance  and  covariances  for  the  involved  processes, 
rY(t)  =  /  eiu;t  RyM  da;, 

rN ,Y(t)  =  Cov(N(0),  Y(t))  =  /  e_ia/t  mw2  {H(a;)-i-l}  RyM  da;, 
rN\Y(t)  =  Cov(N'(0),  Y(t))  =  /  e-iu;t  mia)3  {H(a;)+1}  RyM  dw, 
and 
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(7^  =  /  m“o/*  |H(o;)+l|“  Ry(o;)  da;, 

=  /  m “u/*  |H(a;)  +  l|“  Ryfu)  da;. 

Now,  let  tk,  k  =  1,2,  ...  be  the  times  of  zero  upcrossings  for  the  process  N(t),  i.e. 
the  potential  jump  times,  and  consider  the  movements  of  spring  and  car  during  the  ficti¬ 
tious  flights,  as  defined  by  (54)  and  (55).  The  Slepian  model  for  Y(tk+t)  after  a  zero 
upcrossing  by  N(tk)  consists  of  a  regression  term  on  N(tk)=0,  and  N'(tk)  which,  being  a 
derivative  at  a  level  upcrossing  in  a  stationary  Gaussian  process,  has  a  Rayleigh 
distribution  with  density 

z<rN,  exp(-z“/2<x^,),  z  >  0.  (57) 

Let  (  be  a  Rayleigh  variable  with  density  (57).  The  direct  model  for  Y(tk+t)  is 

then 

Y*(t)  =  b*(t)  +  C0bJ(t)  +  A*(t),  (58) 

where,  with  E(N(0))  =  -gm, 

b*(t)  +  z0bj(t)  =  E(Y(t)  |  N(0)=0,  N'(0)=zo). 

Further,  the  Gaussian  residual  process  A*(t)  is  independent  of  and  has  mean  zero 
and  covariance  function  given  by  the  conditional  covariance  function  of  Y(t)  given  N(0) 
and  N'(0).  However,  if  we  want  to  describe  the  extension  of  the  flight  starting  at  tk  we 
need  also  the  joint  distribution  of  the  derivative  N'(tk),  the  starting  values  in  (51)  and 
(55)  and  the  residual  process  in  (58). 

We  first  introduce  some  notation.  Define  the  random  vectors  v  =  (N(0),  N'(0))  , 

7r  =  (X(0)+Y(0),  X'(0)+Y'(0),  X(0))T  and  note  that  (u,  n)  has  a  5-dimensional  normal 
distribution  with  mean 


=  (-gm,  0,  -gm/k,  0,  -gm/k) 


T 


and  covariance  matrix 
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The  conditional  distribution  of  jr|i/=(0,z)  is  therefore  tri-variate  normal  with  mean 


r  /.  -2  r  -2  , 

— gm/k  +  gmcr  a  a  ,  a  ,  . 

®  ®  n  n,x+y  n  n  ,x+y 
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=  A  +  z  -  B,  say, 


and  covariance  matrix 

E  =  E  -E  E“]E  . 

KK  KV  VV  VK 

To  get  starting  values  for  X,  Y,  X1,  Y'  in  (51)  and  (55),  let  Co  he  a  Rayleigh 

T 

variable,  with  density  (57),  and  let  C  =  (Ci>  (2,  C3)  be  a  3-dimensional  normal  variable, 
independent  of  Co,  with  mean  zero  and  covariance  matrix  E.  Then  the  starting  values 
(X(tk)+Y(tk),  X'(tk)+Y'(tk),  X(tk))T  are  distributed  as 

A  +  CoB  +  C  •  (59) 

The  complete  Slepian  model  for  Y(tk+t)  is  now 


Rjoad(t)  =  b(t)  +  E  Cjbj(t)  +  A(t), 
j=0 

where  the  regression  is  defined  as 


b(t)+  E  Cjbj(t)  =  E{Y(t)  |  v  =  (0,Co)  ,  *  =  A  +  CoB+C}, 
j=0 

and  the  residual  A(t)  is  independent  of  (Co,  Ci,  C2,  C3)  and  has  the  covariance  function 


r^(s,  t)  =  Cov(Y(s),Y(t)  |  v,  k). 
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Equations  (54)  and  (55)  then  give  the  following  expressions  for  the  length  of  the 
spring,  height  of  the  car,  and  height  of  the  wheel  during  a  flight, 

Spring(t)  =  (a3+Cob3+C3)-exp(-kt/c), 

Car(t)  =  PQ  +  (ai  +  (obi+Ci)  +  t(a2+(ob>+0>)  —  gt“/2, 

Wheel(t)  =  Car(t)  —  PQ  —  Spring(t). 

Combining  these  with  the  road  model  (60),  we  get  the  height  of  the  wheel  above  road  as 
Height(t)  =  Wheel(t)  -  Road(t)  = 

=  (ai+Cobi+Ci)  +  t(a2+£ob2+C2)  —  g t~/2  —  (a3+<;*ob3+C3)  exp(— kt/c)  — 

3 

-b(t)  -  E  Cjbj(t)  —  A(t)  = 
j=0 

=  {ai+a2t-gt2/2-a3-exp(-kt/c)-b(t)}  +  <o{bi+b2t-b3-exp(-kt/c)-b0(t)}  + 

+  Ci{  1 — bj(t)}  +  C2{t — b2(t)}  —  C3{exp(-kt/c)+b3(t)}  —  A(t). 

This  model  process  is  of  the  form  (26),  containing  one  deterministic  function,  four  random 
coefficient  functions,  and  one  Gaussian  residual. 

We  now  consider  the  first  time  the  wheel  hits  the  ground  after  a  jump,  i.e.  the  first 
zero  (down)crossing  of  Height(t),  and  we  want  to  find  its  distribution  using  the  method  of 
Section  4.  One  then  has  to  observe  that  Height(t)  already  contains  four  random  coeffi¬ 
cients,  which  have  to  be  conditioned  on.  Further,  to  get  good  accuracy  in  the  approxi¬ 
mation  one  may  have  to  condition  on  some  values  of  the  residual  A(t),  t  hus  increasing  the 
number  of  regressors  to  a  point  where  the  computation  time  becomes  prohibitive. 

It  is  therefore  more  efficient  to  slightly  reformulate  the  model,  including  the  three 
normal  coefficients  Ci>(2>(3  into  the  residual,  and  write  the  model  as  the  sum  of  one  deter- 
minstic  function,  one  Rayleigh  term  and  one  Gaussian  residual, 

Height(t)  =  {ai+a2t-gt2/2-a3-exp(-kt/c)-b(t)}  + 

+  <o{bi+b2t-b3exp(-kt/c)-bo(t)}  +  A°(t),  (61) 

where  now  A°(t)  is  nonstationary  Gaussian  with  a  covariance  function  given  by  r^  and 
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the  covariances  between  0,  (2,  (3-  This  will  permit  the  computer  algorithm  to  pick  as 
regressors  those  values  of  A°(t)  which  contribute  the  most  to  the  distribution  of  the  first 
zero  crossing.  Note  that  (61)  is  of  the  same  simple  form  as  the  original  Slepian  model  (11). 

The  following  numerical  example  shows  the  regression  approximation  for  the  time  of 
first  zero  of  Height(t)  when  the  Gaussian  road  process  Y0(s)  is  of  the  standard  form  with 
spectral  density 

Ry  (a;)  =  <7y/2v/3  for  |u>|  < 

(  .  ais  is  rather  different  from  real  road  spectra  (see  Lindgren  (1981)  for  references)  but  is 
chosen  here  as  a  standard  process.)  Further  k  =  0.64  N/m,  c  =  0.08  Ns/m,  and  m  = 

0.01  kg.  We  have  chosen  two  different  speed  values,  v  =  2  in/s,  4  m/s.  With  try  =  0.5,  the 
jump  intensity,  i.e.  the  mean  number  of  zero  upcrossings  per  time  unit  for  N(t),  is  0.035 
and  0.85,  respectively. 


Figure  5:  Regression  approximations  to  jump  length  density  for  v=2  and  v=4; 
( -  k  =  0,  *  —  *  —  *  k  =  2) 


With  a  model  of  the  form  (61),  the  regression  approximation  of  the  time  for  the  first 
zero  are  shown  in  Figure  5  for  the  two  v-valucs.  For  v  =  2  m/s  it  suffices  to  use  one  extra 
regressor,  while  for  v  =  4  m/s  one  needs  two  in  order  to  find  the  characteristic  details  of 
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the  distribution.  As  seen,  the  speed  is  important;  for  v  =  4  m/s  the  distribution  of  jump 
length  has  two  tops,  while  it  is  unimodal  for  low  speed.  This  can  be  explained  by  the 
periodicity  of  the  road  profile.  For  v  =  4  m/s,  the  wheel  can  jump  over  the  first  hill  and 
hit  the  second. 


9 

5.2  Regression  approximations  for  crossings  in  the  ^ processes 
2 

A  \  -process  is  the  sum  of  squares  of  Gaussian  processes, 

9  n  9 

r(t)=  e  zT(t), 

i=l  1 

where  Z^(t),  ...,  Zfi(t)  are  independent,  with  mean  zero  and  variance  one.  In  the  general 
2 

\  —process,  Zj(s)  and  Z^(t)  need  not  be  independent  fcr  i  ^  k,  s  #  t,  but  in  order  to  keep 

9 

complexity  down,  we  shall  here  consider  only  ^“-processes  with  independent  components 
with  common  covariance  function  r(t).  Let  A9  =  V(Z!(t))  =  -r"(0).  (The  envelope  of  a 

9 

Gaussian  process  is  a  simple  example  of  a  \“-process  with  dependent  components;  see 
Ditlevsen  &  Lindgren  (1988),  and  Lindgren  (1989)  for  details  on  the  Slepian  model.) 

The  Slepian  model  for  \  ( • )  after  upcrossings  of  a  level  u“  was  studied  by 
Aronowich  &  Adler  (1986).  In  the  following  simple  form,  given  by  Lindgren  (1989), 

Y>)  =  (b0(t)-<,b,(t)+A,(t))2  +  £  K^M+AjM)2, 

Cj  is  a  standard  Kaylcigh  variable,  •••»  Cn  are  independent  standard  Gaussian  vari¬ 
ables  and  the  residual  processes  Aj  are  independent  zero  mean  Gaussian  with  covariance 
function  given  by  (8).  Note  that  here  are  ony  two  types  of  b-functions,  given  by 

bQ(t)  =  ur(t), 

bj(t)  =  r'(t)/A^2. 

Suppose  we  are  interested  in  the  length  T  of  an  excursion  above  u  by  the  \“- 

2  4-9 

process  x  (0>  Le.  of  the  first  T  for  which  Y^(T)  =  u  .  In  order  to  put  the  problem  in 
the  general  setting  of  Section  4,  we  introduce  the  process 
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flt)  =  Y+(l)-u2  = 

=  (yO-^bpRAp))2  +  S  U-ib,(t)+Aj(t))2-u2  (62) 

As  in  Section  4,  the  regression  method  is  based  on  a  decomposition  of  £(t)  into  a  re¬ 
gression  term  and  independent  residual  processes,  but  here  the  decomposition  is  not  linear 
as  in  equation  (26)  but  quadratic,  which  leads  to  mixed  terms  of  the  form  (iAj(s). 

However,  the  most  important  common  property  of  (26)  and  (62)  is  that  in  both 
cases,  the  ^-process  is  defined  by  a  mapping  of  (  =  ((p  ...,  Cn)  and  A(t)  =  (Aj(t),  ..., 
An(t)),  i.e.  £(t)  =  F(t,  (,  A(t)).  Under  some  additional  assumptions  on  F,  one  can  gene¬ 
ralise  the  results  of  Section  4,  to  cover  this  more  complicated  situation;  see  Rychlik  (1987b, 
1989b). 

Like  in  Section  4,  the  regression  approximation  of  T  is  defined  as  the  correspon¬ 
ding  variables  in  the  regression  functions  obtained  by  replacing  in  (62),  the  residual 
processes  Aj  by  the  regressions  Aj  on  some  additional  explanatory  vectors  Y*  =  (Yj,..., 
Y‘),  i.e 

{£(0  =  (b0(t)-C,bt(t)+A,(t))2  +  E  (CjbpHAjd))2-^,  (63) 

where 

Ai(t)  =  E(Aj(t)|Yi).  (61) 

In  general,  the  number  of  regressors  in  the  vector  V  may  depend  on  i,  i.e.  k=k*. 
However,  in  our  examples  the  residual  processes  Aj  are  independent  and  identically 
distributed,  and  we  choose  the  same  function  to  define  the  regressors  for  each  residual,  i.e. 
Y1  =  Y(Aj(-)).  This  choice  leads  to  independent  and  identically  distributed  regressions 
(64).  In  addition,  since  the  first  squared  process  in  (62)  is  exactly  the  same  as  the  Slepian 
model  (11),  we  can  use  the  same  procedure  to  select  Y1  =  Y( Aj(  • ))  as  we  used  in  Section 
4.6  for  the  distribution  of  excursion  length  distribution. 

r  r  j* 

We  shall  denote  the  regression  approximation  of  T  by  T^,  i.e.  =  T(^).  Since 
Aj  are  Gaussian  processes,  the  regressors  Y1  can  be  choosen  using  the  same  methods  as  in 
Section  4.5. 
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We  begin  with  the  true  density  of  T.  As  before,  using  the  implicit  definition  of  T, 
we  find  a  transformation  which  represents  one  of  (Vs  as  a  function  of  T,  the  remaining  (Vs 
and  the  residual  processes,  Aj,  i.c. 


>  9 


(b0(T)— CjbjfTj+A^T))-  +jM(ib1(T)+Aj(T))“  -  u“  =  0.  (65) 

Since  ^  has  a  larger  variance  than  the  other  ^-variables,  we  solve  (65)  for  (j.  Assuming 
bj(t)  +  0,  then,  for  T  =  t,  ( C>*--Cn)  =  x  and  A  =  u,  we  have 

Cl  =  Hu2  -  MxjbjfD+^d))2}1/2  -  b0(t)  -  ^(t))/b,(t)  =  q2(t,x,u;)  (66) 


C{  =  ({u2  -  £  (xjb1(t)  +  o;i(t))2}1/2  -  b0(t)  -  ^(OJ/bjd)  =  qi(t,x,o;).  (67 

i — 2 

Similarity  to  (39),  the  Jacobians  Jj,  i=l,2,  of  the  transformations  qjt  (66,  67),  are 
obtained  by  the  implicit  function  theorem  for  evaluation  of  d^/dT,  viz. 


Jj(t,X,Cj)  =  - 


35?lxJs) 


(-l)i(u2-E?=2(xib1(t)+Wi(t))2}1/'2b1(t) 


where  the  function  C[xJs)  is  defined  by  "£(s)|  £(t)=0,  C2,-,Cn=x,  A=u/\  i.e. 


CuJs)  =  (bQ(s)— qj(t,x,o>)bj(s)+u;1(s))“  +  E  ^(xjb^sHu^s))2  -  u2.  (69) 

Then  the  true  conditional  density  of  (T|  £,,...,£  )  can  be  written  as  an  expectation 
over  all  sample  paths  of  A  =  (Alv..,An),  of  the  density  (T|  <0,...,Cn,A),  which  is  defined 
by  the  variable  transformation  (66,  67),  i.e. 

2 

fT|C2-Cn^t,X^  =  -5  EA  lJi(t’x’A)hj(t^)exP{-°-5qi(t,x,A)2}  ,  (70) 

1—1 

since  Cj  is  a  Rayleigh  variable  independent  of  (2,  •••,  Cn- The  indicator  function  I(t;  • )  is 
defined  by  (34),  i.e. 


I(t;f)  = 


1  if  f  ( s)  >0  for  all  s,  0<s<t, 
0  otherwise . 


Once  again,  the  density  of  T  can  be  written  in  the  "marked  crossing"  form,  see  (42),  i.e. 
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fT(t)  =  E(I(t;OU'(t)|  |  «t)=0).f^(t)(0). 

The  expectation  in  (70)  is  difficult  to  evaluate  exactly,  since  the  indicator  I(t;£tx,J 
depends  on  the  whole  realization  of  £[xty,  i.e.  on  A. 

Finally,  the  conditional  density  of  the  regression  approximation  t£|  ^2 » -  *  * » Cn  ob¬ 

tained  by  replacing  in  (70)  the  residual  processes  A  =  (Aj,...,An)  by  their  regressions 
A  =  (A1V..,A  )  as  in  (63).  In  the  simplest  case,  k=0,  the  residuals  Aj  are  replaced  by 
constant  functions  equal  to  zero.  Then  the  expectation  in  (70)  disappears,  leading  to  an 
explicit  density,  viz. 


c,2...U 


( 1 1 x )  =  $  l(t;^[xo)- ! Jj(t,x,0)|qi(t,x,0)exp{-0.5qi(t,x,0)2}j. 


0  00  1.00  2  00  3  00  <  00  0  00  1.00  2.00  3.00 

Figure  6:  Regression  approximation  of  excursion  length  density  after  a  u— upcrossing  for 

o 

\“-process,  the  sum  of  squares  of  two  standard  normal  processes;  (a)  u  =  2,  (b)  u  =  3; 


k  =  0  ( - ),  k  =  1  (  * - *),  k  =  3  ( - ). 
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